﻿/*
 * SphericConformalTransform
 * 
 * initial implementation : 6.7.2011
 *
 * Copyright (c) Adam Sida, asida@centrum.cz.
 * All rights reserved.
 * 
 * for license see Petra.Gis.License.txt
 */
using System;
using System.Collections.Generic;
using System.Runtime.Serialization;
using Petra.Gis.Geometry.Support;


namespace Petra.Gis.ReferenceSystemServices.Transformations
{

  /// <summary>
  /// This is pseudo-datum transform especially used in european coordinate systems
  /// (Swiss, Hungary, Czech, Slovak) where standard cartographic projection (conic, cylindrical)
  /// is at general (oblique) position on spheroid. This projection is one more step 
  /// that must be done in transform between geographic space (lambda, phi) and specific projection (x,y).
  /// Projection uses gauss sphere with radius given by specific standard parallel.
  /// Sphere position relative to Spheriod is also given by oblique standard parallel,....
  /// This projection supports also pole reprojection (cartographic pole). 
  /// New pole position is given by its position in gauss spheric space (lambda, phi).
  /// Transform can threfore produce not only standard (lambda, phi) coordinates, but also
  /// non standard (rotated) spheric coordinates. Coordinate names are defined to U (longitude), V (lattitude).
  /// Projection has different equations for forward and reverse direction.
  /// Forward direction is defined as gauss sphere to ellipsoid.
  /// </summary>
  [DataContract]
  public class SphericConformalTransform1 : InversibleParametricTransform
  {

    /// <summary>
    /// hotine version, temporary
    /// </summary>
    private bool hotine = true;

    //transform params
    private double centralMeridian;


    /// <summary>
    /// Constants used in the transformation.
    /// </summary>
    private double B, A, E;


    /// <summary>
    /// Convenience values equal to A/B, A*B, B/A
    /// </summary>
    private double ArB, AB, BrA;

    
    /// <summary>
    ///<var>v</var> values when the input latitude is a pole. 
    /// </summary>
    private double v_pole_n, v_pole_s;

 
    //transform input parameters
    private double latitudeOfCentre, longitudeOfCentre;
    private double azimuth;

    /// <summary>
    /// computed angle
    /// </summary>
    private double gamma0;
    private double cosgamma0, singamma0;


    private double e;  // eccentricity
    private double e2; // eccentricity squared

    private double u_c;

    //equation variables
    //private double tmpLat;
    //private double sinPhi;
    //private double cosPhi;
    //private double esf;
    
    //oblique constants
    private bool isOblique;
    
    //equation variables
    private double Q, S, V, U, Qp, Sp, Vp, Up, tmp;
    private double u, v;

    AngularUnit srcUnit;
    AngularUnit trgtUnit;




    #region constructors

    /// <summary>
    /// xaml constr.
    /// </summary>
    public SphericConformalTransform1() 
      :base()
    {
      Metadata = new Metadata("SphericConformalTransform", "none", -1, "", "", "provides transformation between Geographics Coordinate systems and gauss sphere");
    }





		public SphericConformalTransform1(GeographicCoordinateSystem targetCS,
                                     double obliqueStandardParallel,
                                     AngularUnit angularParamUnit,
                                     bool isInverse)
			:this(targetCS, angularParamUnit, isInverse,
            new TransformParameters(new XParam("oblique_standard_parallel", obliqueStandardParallel)) )
    {  }





    protected SphericConformalTransform1(GeographicCoordinateSystem targetCS,
                                       AngularUnit angularParamUnit,
                                       bool createInverse,
                                       TransformParameters parameters)
      : base(GeographicCoordinateSystem.CreateGaussSphereCS(targetCS, parameters["oblique_standard_parallel"] * angularParamUnit.RadiansPerUnit),
             targetCS, null, angularParamUnit, createInverse, parameters)
    {
      if (parameters.Count < 1)
        throw new ArgumentException("SphericConformalTransform requires at least 1 parameter");

      if(targetCS.Unit == null)
        throw new ArgumentException("coordinate system must have unit defined");

      Metadata = new Metadata("SphericConformalTransform", "none", -1, "", "", "provides transformation between Geographics Coordinate systems and gauss sphere");
    }


    #endregion constructors




    protected override void updateTransform()
    {

      srcUnit = (AngularUnit)Source.Unit;
      trgtUnit = (AngularUnit)Target.Unit;

      latitudeOfCentre = Parameters["latitude_of_center"] * base.AngularParameterUnit.RadiansPerUnit;
      longitudeOfCentre = Parameters["lontitude_of_center"] * base.AngularParameterUnit.RadiansPerUnit;

      //eccentriticy
      e = ((GeographicCoordinateSystem)Target).Datum.Ellipsoid.Eccentricity;
      e2 = ((GeographicCoordinateSystem)Target).Datum.Ellipsoid.e2;

      double var_1e2 = Math.Sqrt(1.0 - e2);

      double sinphi0 = Math.Sin(latitudeOfCentre);
      double cosphi0 = Math.Cos(latitudeOfCentre);
      double tmp = cosphi0 * cosphi0;

      B = Math.Sqrt(1.0 + e2 * (tmp * tmp) / (1.0 - e2));

      double con = 1.0 - e2 * sinphi0 * sinphi0;
      
      A = B * var_1e2 / con;
      
      double D = B * var_1e2 / (cosphi0 * Math.Sqrt(con));
      double F = D * D - 1.0;

      if (F < 0.0)
        F = 0.0;

      else
      {
        F = Math.Sqrt(F);
        if (latitudeOfCentre < 0.0)
          F = -F; // Taking sign of 'latitudeOfCentre'

      }

      F = F += D;
      E = F * Math.Pow(ProjMath.Tsfnz(e, latitudeOfCentre, sinphi0), B);

      //resolve azimuth
      azimuth = Parameters["azimuth"] * base.AngularParameterUnit.RadiansPerUnit;

      if ((azimuth > -1.5 * Math.PI && azimuth < -0.5 * Math.PI)
          || (azimuth > 0.5 * Math.PI && azimuth < 1.5 * Math.PI))
        throw new ArgumentOutOfRangeException("azimuth out of range");

      gamma0 = Math.Asin(Math.Sin(azimuth) / D);
      // Check for asin(+-1.00000001)
      tmp = 0.5 * (F - 1.0 / F) * Math.Tan(gamma0);
      if (Math.Abs(tmp) > 1.0)
      {
        if (Math.Abs(Math.Abs(tmp) - 1.0) > ProjMath.EPSLN6)
          throw new ArgumentOutOfRangeException("tolerance exception");

        tmp = (tmp > 0) ? 1.0 : -1.0;
      }

      centralMeridian = longitudeOfCentre - Math.Asin(tmp) / B;

      singamma0 = Math.Sin(gamma0);
      cosgamma0 = Math.Cos(gamma0);

      ArB = A / B;
      AB = A * B;
      BrA = B / A;
      
      v_pole_n = ArB
              * Math.Log(Math.Tan(0.5 * (Math.PI / 2.0 - gamma0)));
      v_pole_s = ArB
              * Math.Log(Math.Tan(0.5 * (Math.PI / 2.0 + gamma0)));

      if (hotine)
        u_c = 0.0;

      else
      {
        if (Math.Abs(Math.Abs(azimuth) - Math.PI / 2.0) < ProjMath.EPSLN)
          // LongitudeOfCentre = NaN in twoPoint, but azimuth cannot be 90 here (lat1 != lat2)
          u_c = A * (longitudeOfCentre - centralMeridian);

        else
        {
          tmp = Math.Abs(ArB * Math.Atan2(Math.Sqrt(D * D - 1.0), Math.Cos(azimuth)));

          if (latitudeOfCentre < 0.0)
            tmp = -tmp;

          this.u_c = tmp;
        }

      }


    }





    protected override void transform(ICoordinate coord)
    {

      coord.lon *= srcUnit.RadiansPerUnit;
      coord.lat *= srcUnit.RadiansPerUnit;

      if (IsInverse)
      {
        //from datum to aposphere/gauss sphere...

        //throw new NotImplementedException();
        
        if (Math.Abs(Math.Abs(coord.lat) - Math.PI / 2.0) > ProjMath.EPSLN6)
        {
          Q = E / Math.Pow(ProjMath.Tsfnz(e, coord.lat, Math.Sin(coord.lat)), B);
          tmp = 1.0 / Q;
          S = 0.5 * (Q - tmp);
          V = Math.Sin(B * coord.lon);
          U = (S * singamma0 - V * cosgamma0) / (0.5 * (Q + tmp));

          if (Math.Abs(Math.Abs(U) - 1.0) < ProjMath.EPSLN6)
              throw new ApplicationException();

          v = 0.5 * ArB * Math.Log((1.0 - U) / (1.0 + U));
          tmp = Math.Cos(B * coord.lon);

          if (Math.Abs(tmp) < ProjMath.EPSLN)
              u = AB * coord.lon;

          else
              u = ArB * Math.Atan2((S * cosgamma0 + V * singamma0), tmp);

        }
        else
        {
            v = coord.lat > 0 ? v_pole_n : v_pole_s;
            u = ArB * coord.lat;
        }

        u -= u_c;
        
        coord.lon = u;
        coord.lat = v;

      }
      else 
      {
        //from aposphere/gauss sphere to datum

        //throw new NotImplementedException();

        Qp = Math.Exp(-BrA * coord.lat);
        tmp = 1.0 / Qp;
        Sp = 0.5 * (Qp - tmp);
        Vp = Math.Sin(BrA * coord.lon);
        Up = (Vp * cosgamma0 + Sp * singamma0) / (0.5 * (Qp + tmp));
        if (Math.Abs(Math.Abs(Up) - 1.0) < ProjMath.EPSLN6) 
        {
            coord.lon = 0.0;
            coord.lat = Up < 0.0 ? -Math.PI / 2.0 : Math.PI / 2.0;
        }
        else 
        {
            coord.lat = Math.Pow(E / Math.Sqrt((1.0 + Up) / (1.0 - Up)), 1.0 / B); //calculate t
            coord.lat = ProjMath.Cphi2(coord.lat, e);
            coord.lon = -Math.Atan2((Sp * cosgamma0 - Vp * singamma0), Math.Cos(BrA * coord.lon)) / B;
        }
    
      }

      coord.lon /= trgtUnit.RadiansPerUnit;
      coord.lat /= trgtUnit.RadiansPerUnit;
    
    }





    protected override bool parametersValid()
    {
      return Parameters.ContainsKey("latitude_of_center")
        && Parameters.ContainsKey("lontitude_of_center")
        && Parameters.ContainsKey("azimuth");

    }





    /// <summary>
    /// must override here, spheric conformal doesnt need parameter distance units (only angular)
    /// </summary>
    /// <returns></returns>
    protected override bool graphCompleted()
    {
      return !(this.AngularParameterUnit == null || !parametersValid());
    }






    public override SRSTransform Inverse()
    {
      return new SphericConformalTransform1(
                      (GeographicCoordinateSystem)Target,
                      AngularParameterUnit,
                      !IsInverse,
                      Parameters);
    }





    public override Dimension TransformDimension
    {
      get { return Dimension.dim2D; }
    }





    ///// <summary>
    ///// creates new instance of Spheric conformal transform from geogr. coordinate system
    ///// </summary>
    ///// <param name="gcs">source georg. coordinate system</param>
    ///// <param name="angParamUnit">transformation parameter angular unit (distance is not relevant for SphericConformal)</param>
    ///// <param name="parameters">transform parameters</param>
    ///// <returns>new instance of SphericConformalTransform or null or throws exception when invalid input parameters</returns>
    //public static SphericConformalTransform1 Create(GeographicCoordinateSystem gcs, AngularUnit angParamUnit, bool createInverse, TransformParameters parameters)
    //{
    //  if (parameters["oblique_pole_latitude"] * angParamUnit.RadiansPerUnit > ProjMath.HALF_PI)
    //    throw new ArgumentOutOfRangeException("latitude must be in 0 - 90 deg range");

    //  return new SphericConformalTransform1(gcs, angParamUnit, createInverse, parameters);
    //}


  }

}
